function y=irf(sysm,e,horizon,cumidx)
% y = IRF(SYSM,E,HORIZON)
%  calculates impulse response
%
% INPUT:
%   E       : normalized initial shock
%   HORIZON : length of horizon for responses
%

% Sungbae An
% Updated: 09/14/2005

% check input
if nargin==3
	cumidx=0;
end

% dimensions
nz = size(sysm.R,1);
[ny,ns] = size(sysm.ZZ);
nx = ns - nz;

% initialize output
y = zeros(ny,horizon+1);

% initial response
z = sysm.R*e;
x = sysm.X.lin*[zeros(nx,1);z];

y(:,1) = sysm.ZZ*[x;z];

% horizon loop
for n = 2:(horizon+1)

	z = sysm.T*z;
	x = sysm.X.lin*[x;z];

	y(:,n) = sysm.ZZ*[x;z];

end		% n loop

if sum(cumidx)>0
	% cumulative response
	y(cumidx,:) = cumsum(y(cumidx,:),2);
end

